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Abstract 

The structure theorem of Hadamard-Zolesio states that the derivative of a shape functional is a distri¬ 
bution on the boundary of the domain depending only on the normal perturbations of a smooth enough 
boundary. Actually the domain representation, also known as distributed shape derivative, is more general 
than the boundary expression as it is well-defined for shapes having a lower regularity. It is customary 
in the shape optimization literature to assume regularity of the domains and use the boundary expression 
of the shape derivative for numerical algorithms. In this paper we describe several advantages of the dis¬ 
tributed shape derivative in terms of generality, easiness of computation and numerical implementation. 

We identify a tensor representation of the distributed shape derivative, study its properties and show how 
it allows to recover the boundary expression directly. We use a novel Lagrangian approach, which is ap¬ 
plicable to a large class of shape optimization problems, to compute the distributed shape derivative. We 
also apply the technique to retrieve the distributed shape derivative for electrical impedance tomography. 
Finally we explain how to adapt the level set method to the distributed shape derivative framework and 
present numerical results. 
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Introduction 

In his research on elastic plates [19] in 1907, Hadamard showed how to obtain the derivative of a shape 
functional J{Q) by considering normal perturbations of the boundary dQ of a smooth set U. This fundamental 
result of shape optimization was made rigorous later by Zolesio [13] in the so-called “structure theorem”. 
When J(U) and the domain are smooth enough, one may also write the shape derivative as an integral over 
dQ, which is the canonical form in the shape optimization literature. 

However, when Q is less regular, the shape derivative can often be written as a domain integral even when 
the boundary expression is not available. The domain expression also known as distributed shape deriva¬ 
tive has been generally ignored in the shape optimization literature for several reasons: firstly the boundary 
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representation provides a straightforward way of determining an explicit descent direction since it depends 
linearly on the boundary perturbation 9 and not on its gradient, secondly this descent direction only needs to 
be defined on the boundary. When considering the domain expression, these two advantages disappear as the 
shape derivative is defined on 12 and depends on the gradient of 9, so that a partial differential equation needs 
to be solved to obtain a descent direction 9 on 12. 

It seems that these drawbacks would definitely rule out the distributed shape derivative, however they turn 
out to be less dramatic than expected in many situations and the domain formulation has other less foreseeable 
advantages over the boundary representation. In this paper we advocate for the use of the distributed shape 
derivative and discuss the advantages of this formulation. 

The boundary representation has the following drawbacks. First of all if the data is not smooth enough 
the integral representation does not exist so that the more general domain representation is the only rigorous 
alternative. Even when the boundary representation exists and has the form g 9 ■ n, it is usually not 
legitimate to choose 9 ■ n = —g on c212 for a descent direction if g is not smooth enough, for instance if 
g G L^((912). Therefore, a smoother 9 must be chosen, which requires to solve a partial differential equation 
on the boundary 012. When taking 9 ■ n = —g is legitimate, it might still not be desirable as this may yield a 
9 with low regularity, in which case one needs to regularize 9 on the boundary as well. In these cases the first 
advantage of the boundary representation disappears. The second advantage of the boundary representation is 
that the perturbation field only needs to be defined on the boundary instead of on the whole domain, reducing 
the cost of the computation. Actually, the distributed shape derivative also has its support on the boundary, 
and may be computed in a small neighborhood of the boundary so that the additional cost is minimal. In 
addition, in most shape optimization applications, g is the restriction of a function defined in a neighborhood 
of the boundary and not a quantity depending only on the boundary such as the curvature. Therefore from 
a practical point of view, g must be evaluated in a neighborhood of 012 anyway. Also, in many numerical 
applications, 9 must be extended to a neighborhood of T or even to the entire domain 12. This is the case for 
level set methods for instance, where the level set function must be updated on 12, or when one wishes to 
update the mesh along with the domain update, to avoid re-meshing the new domain. The distributed shape 
derivative then directly gives an extension of 9 well-suited to the optimization problem. 

Recent results have shown that the distributed shape derivative is also more accurate than the boundary 
representation from a numerical point of view; see [27] for a comparison. Indeed functions such as gradients 
of the state and adjoint state appearing in the distributed shape derivative only need to be defined at grid points 
and not on the interface. Therefore one avoids interpolation of these irregular terms. This is particularly useful 
for transmission problems where the boundary representation requires to compute the jump of a function over 
the interface, a delicate and error-prone operation from the numerical point of view. 

Having considered these equivalent expressions of the shape derivative (i.e. boundary and domain expres¬ 
sion) leads to a general form of the shape derivative using tensors. We introduce such a tensor representation 
in Section 3 which covers a large class of problems and in particular contains the boundary and domain 
expression. We show how this abstract form allows to identify simple relations between the domain and 
boundary expressions of the shape derivative. 

In this paper we also extend and simplify the averaged adjoint method from [42], a Lagrangian-type 
method which is well-suited to compute the shape derivative of a cost function in an efficient way. Lagrangian 
methods are commonly used in shape optimization and have the advantage of providing the shape derivative 
without the need to compute the material derivative of the state; see [4,9, 12,37,42,43]. Compared to these 
known shape-Lagrangian methods, the averaged adjoint method is fairly general due to minimal required 
conditions. The assumptions are for instance less restrictive than those required for the theorem of Correa- 
Seeger [12], therefore it can be applied to more general situations such as non-convex functionals. As the 
direct approach our method can also be applied for problems depending on nonlinear partial differential 
equations. In this paper we give an example of application to a transmission problem (in electrical impedance 
tomography - see Section 5). Our method provides the domain expression of the shape derivative and the 
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boundary expression can be computed easily from the tensor representation of the domain expression. 

To complete the numerical implementation aspect, we also show how the domain expression of the shape 
derivative can be used in the level set method framework [3, 15-18,23,24, 35]. The level set method can be 
modified to use the domain expression which leads to a method which is actually easier to implement. Com¬ 
bining all these techniques, we obtain a straightforward and general way of solving the shape optimization 
problem, from the rigorous theoretical computation of the shape derivative to the numerical implementation. 

In Section 1 we recall the concept of shape derivative and the structure theorem on an abstract level. In 
Section 2 a shape-Lagrangian method, the averaged adjoint method, is described. In Section 3 we identify a 
general tensor representation of the shape derivative, establish some of its properties, and give a few examples. 
In Section 4 we explain how to compute descent directions for the distributed shape derivative for use in 
gradient methods. In Section 5 we apply the results of Sections 2 and 3 to the inverse problem of electrical 
impedance tomography. In Section 6 we extend the level set method to the case of the distributed shape 
derivative and hnally in Section 7 we show numerical results for various problems including the problem of 
electrical impedance tomography. 

1 The structure theorem revisited 

Our aim in this section is to describe properties of the shape derivative on an abstract level and to emphasize 
that all representations of the shape derivative satisfy the same structure theorem. 

Let V{D) be the set of subsets of 77 C R'^ compactly contained in D, where the so-called “universe” 
D C R'^ is assumed to be open and bounded. Dehne for /c > 0 and 0 < a < 1, 

C'c^’“(L>,R'^) := {(9 G C'=’"(D,R‘^)| 0 has compact support in D}. (1.1) 

Also for given domain Q C D with at least a boundary we introduce the space of vector held 

:= {0 G C^"(77,R'^)| 0 • n = 0 on 512} (1.2) 

where n is the outward unit normal vector to Q. 

Consider a vector held 0 G Cc'^{D, R*^) and the associated how : 77 —R^, t G [0, r] dehned for 
each xo G 77 as 4*^ (xq) := x(f), where x : [0, r] —^ R solves 

x{t) = 6{x{t)) fortG(0,r), x(0) = xq. (1.3) 

We will sometimes use the simpler notation when no confusion is possible. Since 6 G Cc’^{D, R'^) 

we have by Nagumo’s theorem [32] that for hxed t G [0, r] the how is a homeomorphism from D into 
itself and maps boundary onto boundary and interior onto interior. Further, we consider the family 

Qt := 4>f(L!) (1.4) 

of perturbed domains. 

In the following let J : ^ R be a shape function dehned on some admissible set ^ C 'P{D). 

Definition 1.1. The Eulerian semiderivative of J at 17 in direction 9 G Cc’^ (77, R'^), when the limit exists, is 
dehned by 

(1,5) 

(i) J is said to be shape differentiable at 17 if it has a Eulerian semiderivative at 17 for all 
and the mapping 

dJ(17) : C^iD, R'^) ^ R, 0 ^ dJ(17)(0) 
is linear and continuous, in which case dJ(17)(0) is called the shape derivative at 17. 
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(ii) The shape derivative dJ{^) is of finite order if there is an integer I > 0 and a constant c > 0 such that 
for each compact K C D 

\dj{n){6)\ < c\\d\\i yeec^{K,-R^), 

where ||0||; := X]|q,|<; |T>"0|oo- The smallest such integer f > 0 is called order of dJ{Q,). 

The shape derivative from Definition 1.1 has a particular structure. Intuitively, it is clear that the form 
functional stays constant for a transformation $ that leaves Q unchanged, that is 4’(D) = D, even if some 
points inside D move and consequently the shape derivative is zero in this case. This property is valid when 
Q is open or closed; cf. [13]. Mathematically, this is expressed in the following basic theorem proved in [44]. 

Theorem 1.2. Let D G ^ be open or closed. Let 6 G Cc'^{D, R*^) be a vector field wifh compacf supporf in 
Q and denote by ifs flow defined in (1.3). Then we have 

dj{n){e) = 0 . 

Nofe fhaf fhe shape derivafive of J{Q) always exisfs for vecfor fields wifh compacf supporf in 0, even if if 
does nof exisf for ofher vecfor fields. An imporfanf consequence of Theorem 1.2, also for numerical mefhods, 
is fhaf independently of the representation of fhe shape derivafive and fhe regularify of fhe domain D, fhe 
values of 9 outside the boundary of D have no influence on fhe shape derivafive. 

Corollary 1.3. Lef D G ^ be a sef wifh C^-boundary. Assume fhaf J is shape differenliable on fjd. Lef 
9 G {D, R'^). Then we have 

dJ{n){9) = 0. 

The previous discussion immediately yields fhe following fundamenfal resulf of shape opfimizafion. 

Theorem 1.4 (Sfrucfure Theorem). Assume T := dO. is compacf and J is shape differentiable. Denofe fhe 
shape derivative by 

dJ(D) : C“(L>,R'^) ^ R, 9^dJ{^){9). (1.6) 

Assuming dJ{Vl) is of order k > 0 and T of class fhen fhere exisfs a linear and continuous functional 
: C'^(r) —> R such fhat 

dJ{n){9) = g{9\r-n), (1.7) 

Proof. See [13, pp. 480-481]. ■ 

2 Shape derivatives via averaged adjoint method 

Lagrangian mefhods in shape opfimizafion allow fo compufe fhe shape derivafive of funcfions depending on 
fhe solution of partial differential equations without the need to compute the material derivative of the partial 
differential equations; see [13] for a description of such a method in the linear case. Here we extend and 
simplify the averaged adjoint method, a Lagrangian-type method introduced in [42]. With this approach the 
computation of the domain representation of the shape derivative is fast, the retrieval of the boundary form is 
convenient, and no saddle point assumptions is required unlike in [13]. 

Let two vector spaces E = E{Q),F = F{Q) and r > 0 be given, and consider a parameterization 
Qt = ^t{^) for t G [0, r]. Ultimately, our goal is to differentiate shape functions of the type J(Ui) which 
can be written using a Lagrangian as J(Dt) = where tt* G E{Q,t) and G F{Q,t)- The 

main appeal of the Lagrangian is that we actually only need to compute the derivative with respect to t of 

C{Qt, f>i f’) to compute the derivative of indeed this is the main result of Theorem 2.1, but this requires 

a few explanations. 
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Since C{Qt, 0, V')^ is often constituted of integrals on $^(0), using a change of variable we can rewrite 
these integrals to integrals on the fixed domain fl, and consequently transfer the dependence on t to the 
integrand. However, in the process appear the composed functions G E{Q) and ^ whose 

derivatives are not straightforward to compute since ip and ip are defined on fhe moving spaces E{Q,t) and 
E{nt). 

Forfunafely, and fhis is fhe crucial poinf of fhe shape-Lagrangian approach, fo compufe fhe shape deriva- 
five we can reparameferize fhe problem by considering 'I't o ip) insfead of £(0^, (p, ip), where 

'I't is an appropriafe bijecfion befween E{Q) and E{Qt), and p G E{^), ip G E{Q,). Now fhe change 
of variable in fhe infegrals yields funcfions p and ip in fhe infegrand, which are defined on fixed spaces. 
In fhis paper E and E are i?^-spaces, and in fhis case we may consider fhe parficular reparameferizafion 

spaces such as Ff(curl; H), ofher fransformafions 'I't can be used; see [21,26,28]. 
Thus we are led fo consider general funcfions of fhe fype G : [0,r] x£xF—)-R wifh 

G{t,p,ip) := /:{<^ti^),po^-^,iPo<^-^). 

This is precisely whaf we do in (5.29) when showcasing an applicafion of fhe mefhod. The main resulf of fhis 
secfion. Theorem 2.1, shows fhaf fo obfain fhe shape derivafive of C, if is enough fo compufe fhe derivative 
wifh respecf fo f of G while assigning fhe proper values fo p and ip. The main ingredienf is fhe infroducfion 
of fhe averaged adjoinf equation. 

In addifion, in fhis paper we consider fhe following specific form 

G{t,p,ip) := a{t,p,ip)+ b{t,p), (2.1) 

where 

a : [0, r] X £ X F —)• R, 6 : [0, r] x F —)• R, 

are funcfions such fhaf ip a{t, p, ip) is linear for all t G [0, r] and p ^ E. The function G is commonly 

called Lagrangian, hence fhe name of fhe mefhod. In fhe applicafions we have in mind, fhe funcfion b arises 

from fhe objecfive funcfion while a corresponds fo fhe consfrainf, affer fransporfing back fo fhe fixed domain 

0 . 

Throughouf fhe paper, fhe Greek leffers p and ip are used for variables, while fhe roman leffers u, p are 
used for fhe solufions of fhe sfafe and adjoinf sfafes, respecfively. 

Lef us assume fhaf for each t G [0, r] fhe equation 

d^G{t, u^, 0; Ip) = a{t, u^, ip) = 0 for all ip e E. (2.2) 

admifs a unique solufion G F. Furfher we make fhe following assumpfions for G. 

Assumption (HO). For every {t,ip) G [0,r] x F 

(i) [0,1] 9 s I—)• G{t, svP + s{v!' — vP),ip) is absolutely continuous. 

(ii) [0,1] 9 s I—)■ d^G{t, su^ + (1 — s)u°, ip; p) belongs to L^{0, l)for all p £ E. 

When Assumption (HO) is satisfied, for f G [0, r] we infroduce fhe averaged adjoint equation associafed 
wifh and Find £ E such fhaf 

dipG{t,su* + {1 — s)u^,p*;p) ds = 0 for all (^gF. (2.3) 

Notice fhaf, in view Assumption (HO), for all t £ [0, r], 

G{t,u^,p‘') — G{t,u^,p^) = [ d^pG{t,sv!'+ {1 — s)u^,p^;v!'— u^) ds = Q. (2.4) 

Jo 

We can now sfafe fhe main resulf of fhis secfion. 
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Assumption (HI). We assume that 

lim ~ 

t\o t 


= dtG{0,u^,p^). 


Theorem 2.1. Let (HO) and (HI) be satisfied and assume there exists a unique solution of the averaged 
adjoint equation (2.3). Then for G F we obtain 


^b{t,u%=o = ^{G{t,u\il!))\t=o = 9iG(0,n°,p°). 


(2.5) 


Proof. Put g{f) := G(f, n*, 0) — G(0, vP, 0), and note that g{t) = G{t, u^, f) — G(0, vP ^ f) for all ■0 £ F 
and g{0) = 0. We have to show that 

Git,up0)-Gi0,uP0) 

g (0) := lim- exists. 

^ ^ t\o t 

Thanks to Assumption (HO) we ean define the averaged adjoint p^ and using that G is affine with respeet to 
the third argument, we obtain 


g{t) = G{t,u ,p)-G{t,u ,p)+G{t,u ,p)-G{0,u ,p). 

' -V-' 

=0 in view of (2.4) 

Dividing by f > 0 and using Assumption (HI) yields 

'rn\ r 9{t)-g{0) ,. G(f,«°,p*) - G( 0 ,u°,p*) 0 On 

»(“> = — 1 — = ,\“o- 1 -= a<G(0.« . p ) 

whieh eoneludes the proof. 

Remark 2.2. In terms of a and b, equation (2.3) reads: 

/ d^pa{t, sv!' + (1 — s)u^,p^-, (p) ds = — / di^b{t, su^ + (1 — s)u^‘, p) ds 

Jo Jo 

for all (0 G F. If 99 I—)■ a{t, p, 0) is in addition linear, then (2.3) beeomes 

a{t, p,p^) = — dipbit, svJ + (1 — s)vP-, p) ds 

Jo 

for all p ^ E. 


3 Tensor representation of the shape derivative 

In this seetion we identify tensor representations of the shape derivative that eorrespond to a large elass of 
problems studied in the literature for PDE-eonstrained shape optimization. This tensor representation has 
several interesting properties that we investigate. In partieular we exhibit the link between this tensor repre¬ 
sentation and the usual boundary expression of the shape derivative. 
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3.1 Definition and properties 

Definition 3.1. Let ri G ^ be a set with -boundary, A; > 1. A shape differentiable function J of or¬ 
der k is said to admit a tensor representation if there exist tensors S; G R'^)) and Gi G 

L^(9n;£^(R'^,R'^)), I = 0,fc, such that 

k „ 

dj{n){e) = V / Si-D^edx + 

1=0 

where D^O := DO — (DOn) (8) n is the tangential derivative of 6 along dQ. Here £*(R'^,R‘^) denotes the 
space of multilinear maps from R'^ x • • • x R'^ to R*^. 

Most if not all examples involving PDEs from classical textbooks [13,20,41] can be written in the form 
(3.1). 


Ian 


Gi ■ 0 ds for all 0 G (D, R' 


(3.1) 


Remark 3.2. (a) A particular case of the tensor representation (3.1) is the Eshelby energy momentum 

tensor in continuum mechanics introduced in [14]; see also [33]. In this case only Si is not zero. 

(b) When J is is shape differentiable in Q. then by definition 9 i—;■ dJ{Q){6) is a distribution, and if OH is 
compact, the distribution 6 i—)• dJ{Q){6) is of finite order. 

(c) If dJ(H) is of order A: = 1 and \dJ{Q){6)\ < C'||0||Hi(t?,R‘*) ^ ^ C^{D,'R^) then by density 

of C^{D, R*^) in Hq{D, R'^) the derivative dJ{Q) extends to a continuous functional on Hq{D, R*^), 
that is, 

|dJ(H)(0)| < for all 0 G «-'")• 

Therefore by the theorem of Riesz, we obtain a vector field W in Hq{D, R*^) such fhaf 

V0 G Ro(L>,R'^), dJin){e)= [ DW ■ DO + W -Odx 

Jd 

and fhis defines a fensor represenfafion wifh Si = DW, So = W, Si = 0 and So = 0. 

(d) The assumpfion fhaf H be a sef of class can be reduced if S; = 0 for all 0 < A:o < ( < A;. 


The fensor represenfafion (3.1) is nol unique in fhe sense fhaf fhere mighf be several ways fo choose fhe 
fensors Si and S;. This is expressed by fhe facl fhaf fhese tensors are correlafed. We describe fhese relafions 
below in fhe case A; = 1 in Proposifion 3.3, which also describes fhe link befween fhe fensor represenfafion 
and fhe usual boundary represenfafion (1.7) of fhe shape derivafive. 

Proposition 3.3. Eef H be a subsef of D wifh C^-boundary. Suppose fhaf fhe derivafive (iJ(H) has fhe 
represenfafion 

dJin){e)= [ Si-DO + So-Odx + [ Gi-D rO + Go-Ods. (3.2) 

Jd Jan 

If Si is of class in H and D \ H fhen indicafing by -|- and — fhe resfricfions of fhe fensors fo H and 
D\Q, respecfively, we gef 


— div(S]*’) -h S[j' = 0 in H 

— div(S]") -h Sq = 0 in 77 \ H. 


(3.3) 


Moreover, we can rewrite fhe fensor represenfafion as a disfribufion on fhe boundary: 


dJ{n){e) = [ [(S+ - S^)n] -e + Gi-DrO + Go-Ods 
Jan 


where n denotes the outward unit normal veetor to 0,. 

If the boundary is and ©i G R^)), then we obtain a more regular distribution, 

the so-ealled boundary expression of the shape derivative: 

dJ{Q){9)= j gid-nds, (3.4) 

Jdn 

where 

gi := [(S)'' — S^)n] • n + ©o • + ®i • D^n — divr(©fn) + T-L{&[n ■ n). (3.5) 

and % = divr n denotes the mean eurvature^ of dO. while divp := tr(L)r) is the tangential divergenee. 

Proof. Applying Theorem 1.2 we have 

dJ{n){e)= [ Si-De + So-9dx+ [ Gi- Dr9 + eo-9ds = 0 for all 0 G Cl(S2 U (T» \ H), R"*). 
Jd Jdn 

An integration by parts shows (3.3). 

Then, when dQ is C^, replaeing (3.3) in the expression of the shape derivative and using Green’s formula 
we obtain 


dJ{n){9) = [ ©1 • Dr9 + ©0 • 0 ds + [ [{Sf - S^)n] • 9 ds 
Jdn Jdn 

+ j" {-div{Sf) + Sq) ■ 9 dx + _(-div(S^) + Sq ) • ddx ^3 

Jo. V J V 

=0 =0 

Wf [{Sf-S^)n]-9 + ei-Dr9 + eo-9 ds for all 9 £ 

Jdn 

With a slight abuse of notation we keep the same notation n for the extension of the normal to a neighborhood 
of dQ. Let 9 G C^{D, R*^) and define 9r ■= 9 — {9 ■ n)n the tangential part of 9. Then 0,- • n = 0 on cAf and 
henee from the strueture theorem we get dJ{Q){9r) = 0 whieh yields in view of (3.6): 

dJ{n){9) = dJ{fl){{9 ■ n)n) 

= f ((S)'” — Sj“)n • n){9 • n) + ©i • Dr{n{9 ■ n)) + (©o ■ n){9 ■ n) ds 
Jdn 

= f {{Sf — S^)n ■ n){9 ■ n) + {&o ■ n){9 ■ n) ds 
Jdn 

+ / Gi • Drn{9 • n) + n ■ 6iVri9 • n) ds, 

Jdn 

where we used that for all funetions / G C'^(R'^, R'^) and g G C'^(R'^) we have 

D{gf) =gDf + f® Vg. (3.8) 

Finally using ©i G £^(R'^, R'^)) we integrate by parts on the boundary dO. to transform the last 

term in (3.7) 

f n • ©iVr(d • n) ds = f (—divr(©fn) + 7((©fn • n))(d • n)ds. 

Jdn Jdn 


(3.7) 


*We define the mean curvature as the sum of the principal curvatures Ki, that is, % := 
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Therefore (3.7) reads 


dJ{i}){6) = j ((S)*" — )n • n){9 ■ n) + (So • n){6 • n) ds 

Jdn 


+ / &i ■ Drn{0 ■ n) + {—divri&in)+'H{&in ■ n)){9 ■ n) ds, 
Jdn 


whieh we ean rewrite as (3.4). 


(3.9) 


Remark 3.4. In Proposition 3.3, if Si = 0, one ean still obtain (3.3) when is only Lipsehitz instead of C^. 

Corollary 3.5. Let the assumptions of Proposition 3.3 be satisfied. Suppose that the tensor Si : OH —)■ 
£(R'^, R*^) has the form Si = a{I — n®n), where a G (^^(clll). Then (3.4) simplifies fo 


dJ{Q){0)= f gi9-nds, (3.10) 

Jdn 

where gi is given by 

gi := [(S)^ — S)“)n] • n + So • n + aH. 

Proof. Firsf & f = a{I — n® n)^ = Si and Sfn = a{I — n® n)n = a(n — (n • n)n) = 0, fhus fhe two 
last terms in (3.5) vanish. Coneerning the third term in (3.5) we write 

Si • Dj-n = a{I — n®n) ■ Dyu = a(tr(i9rn.) — {n®n) ■ Hpn) = a(divrn — {Dyuti) ■ n) = oH. 

where we have used {DYun) ■ n = 0. ■ 


Remark 3.6. The partieular tensor Si = a(/ — n (8) n) in Corollary 3.5 is eommonly eneountered in shape 
optimization problems. In faet, (3.10) eorresponds to a standard formula that ean be found in most textbooks 
on shape optimization. 

Remark 3.7. Reeall that for given veetor fields 9, C, the seeond order shape derivative is defined by 

d^jimm := j^dj{^Hm0)\t=o. 

Onee we have identified a tensor representation (3.1) for the shape derivative dJ{Q){9) for fixed 9, if is 
eonvenient to differentiate it to also obtain a tensor representation for the seeond-order shape derivative. 
Further, Proposition 3.3 or Corollary 3.5 ean also be applied to obtain a boundary expression for the seeond 
order shape derivative. 


Similar relations as in Proposition 3.3 eould be obtained for any tensor representation of order k. For 
instanee in the ease /c = 2 we obtain the relations 


^S+-div(S+) + S+= 0 inn, 
. 4 S 2 — div(S)") + Sg = 0 in D\n, 


(3.11) 


where (. 482 ); = Ei^,j=i dxiXj{S 2 )iji- 

Using the averaged adjoint approaeh from Theorem 2.1 yields the tensor representation (3.1) of the shape 
derivative. Then Proposition 3.3 ean be used to immediately derive the standard boundary expression of the 
shape gradient from this tensor representation. 
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3.2 Examples of tensor representations 

In this section we present several examples of representations corresponding to Definition 3.1 and apply the 
observations from Section 3.1. 

First order tensor representation 

A basic example of a first order tensor representation of the shape derivative is for 

J(D) = f f dx + f gds 

Jn Jan 

with f,g£ C'^(R'^). Then one easily computes 

dJ{^}){9)= f Vf ■ 6 + f div{9) dx + f Vg ■ 9 + g div-p 9 ds. 

Jn Jan 

The corresponding tensor representation (3.1) is 

S+:=//, Sr:=0, S+:=V/, Sq := 0, &i := g{I - n ^ n), 6o := V 5 . 

Note that ©1 has the form assumed in Corollary 3.5. Applying this Corollary, assuming the domain has 
enough regularity, we obtain in view of (3.10) the classical formula: 

dJ(n)(0) = f gi9-nds, 

Jan 


where gi is given by 


Note that in the particular case 


91 ■= f + dng + g'H. 

/ = 0 we have obtained as a byproduct the formula 
Vg ■ 9 + g divp 9 ds = / {dng + gJi) 9 ■ n ds, 


(3.12) 


Jan Jan 

and when in addition dng = 0 or is defined only on OD, (3.12) becomes fhe classical tangential Green’s 
formula-, see for instance [20, proposition 5.4.9]. 


Non-homogeneous Dirichlet problem 

The following problem was already considered for instance in [13]. Here we present a fairly easy way to 
compute the shape derivative. Let D be an open and bounded subset of that is contained in an open and 
bounded set D. Consider 


-Av = fmn, (3.13) 

V = g on on, (3.14) 

where / G L‘^{D) and g G H‘^{D). Introducing the variable u := v — g, we observe that (3.13)-(3.14) is 
equivalent to the homogeneous Dirichlet problem 


Au = Ag + / in D, 
= 0 on dQ. 


(3.15) 

(3.16) 
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Consider the cost function 


J(n)= / \v-Ud\‘^dx= / \u + g-Ud\‘^ 
Jn Jn 

The weak formulation of (3.15),(3.16) reads: 


dx. 


(3.17) 


Find ?x G : [ Vu-V^ljdx= [ -Vg ■ V'ljj + f'ljj dx for all^P £ (3. 

J Q J Q 


18) 


Note that the previous weak formulation is already well-defined for arbitrary open and bounded set Q. We do 
not need to impose any regularity on Q. The perturbed problem of the previous equation, which is obtained 
by considering (3.18) on <^^((7) and performing a change of variables, reads: find such that 

[ A{t)Vu^-VtPdx = [-A{t)Vg* ■Vij + ^{t)f'tPdx for all V'G (3.19) 

J Q J Q. 

where ^{t) := det(Z7<h4) and A{t) := . The following continuity result is standard: 

Lemma 3.8. There exists a constant c > 0 such that \\u^ — ^ ^ [Oj t]. 

Introduce 

:= / A{t)V(f ■ Vtp dx + / A{t)Vg^ dx 

Jn J n 

b{t,(p) := [ ^{t)\<f + g^-u^Xdx. 

Jn 

where g^ := g o and := Ud o 

Recall that the associated Lagrangian (2.1) is G{t,(p,'ip) = a{t,ip,ip) + b{t,ip). The averaged adjoint 
equation (2.3) reads 

/ A{t)Vip ■ Vdx = / {vP + vP + 2 g*'— 2 v!'d)ip dx for all y? G T7d(^)- 
Jn Jn 

The following continuity result for the adjoint is standard: 

Lemma 3.9. There exists a constant c > 0 such that 

\\P^ - P^\\m{n) < ct for allf G [0,r]. 

One readily verifies that all assumptions of Theorem 2.1 are satisfied, except for (HI) which we now 
prove. Indeed using p* —)■ p^ in 77d(^) ^ 1° ll*® strong differentiability of f i—)• A{t) and 

1 1 -). we get 


G{t, u^,p^) — G{0, u^,p^) 


lim 

t\o 


0 


= lim 

+ 




- I 


Vtt° • Vp* dx + 


\Jn\ t / Jn 

f -\u^ + 9- Ud? 


P A{t)Vg^-Vg 


Vp* dx — 




m 


= dtG{f)yy), 


which shows that (HI) is satisfied. 
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Hence, applying Theorem 2.1 yields 

dJ(H)(6*) = dta{0,u,p) + dtb{0,u), 


which is by definition equivalent to 


dJ{n){e)= [ A'{0){Vu-Vp + Vg-Vp)dx+ [ V(Vp • 0) • Vp - div(/0)pda 
Jq Jo. 

+ / d\Y{9)\u +g - Ud\^ + {y{g - Ud) ■ 0)){u +g - Ud)dx. 

Jn 

Since A'(0) = (div 6)1 — DO"'" — D9 we obtain the tensor representation (3.2) with: 


(3.20) 


Si = I{Vu ■ Vp + Vp • Vp — fp + \u + g — Udl"^) — Vm 0 Vp — Vp 0 Vtt — Vp 0 Vp, 
So = D'^gVp-pVf + {u + g - Ud)V{g - Ud), 

©1 = 0 , ©0 = 0 . 


Now applying (3.10) we get immediately 

= Vm • Vp + Vp • Vp - /p + |rt + P - ttrfp - 2dnUdnP - dnpdng. 

Using the definition of the tangential gradient and p = 0, u = 0 on T implies Vrtt = Vrp = 0, so we obtain 
the simpler expression 


91 = -dnudnP +\u + g - Ud\‘^ = -dn{v - g)dnP + Iv - Mdp. 

Finally, substituting back u = v — g we obtain the formula 

dJ{D){9)= [ {-dniv - g)dnP+\v - Ud\‘^)9 ■ nds. 

Jan 

This formula can be found for instance in [13, p. 566, Formula 6.38], where the adjoint has the sign opposite 
to our case. 

Elliptic problem: first order tensor representation 

Suppose that Q C D C R'^ is a smooth bounded domain, where D C is the smooth “universe”. Let us 
consider the Dirichlet problem: 


— div(MV?x) + ti = / in U, 
ti = 0 on (9U, 

where M G is a positive definite matrix. Consider the cost function 

J(U) = / \u — Ud\^dx^ 

Jn 

where Ud G C^(R‘^). Let us introduce 

:= [ {MQ'Vif ■ Q'Vijj + dx — j ^{t)f'ijjdx 

Jn Jn 

b{t,<p) := [ C{t)\p - u'dl"^ dx, 

Jn 


(3.21) 


(3.22) 
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where Q* := ^ and ^(t) ;= det(D<I>t). Then the weak formulation of (3.21) on the perturbed domain 

Qt, onee transported baek to is 


a{t,u^,'ifj) = 0 for all-0 S-f^o (^^)- 

The Lagrangian eorresponding to the minimization of J{Q) and the PDE eonstraint (3.21) is 

G{t,ip,Tp) = b{t,(p) + (3.23) 

It ean be shown using Theorem 2.1 that dJ{^){6) = dtG{0,u,p), where p G H^i^) denotes the adjoint 
state: 


MVijj ■ Vp + pip dx 



Ud)ip dx 


for all Ip G i^Q (12). 


(3.24) 


The tensor representation (3.1) of the shape derivative of J( 12) in direetion 9 G G^{D, R'^) is given by 


dj{n){e)= [ Si- De + So-edx. ( 3 . 25 ) 

Jn 

where we use the relation (Vp (8> Vn) • DO = DOVu ■ Vp to get the tensors 

50 = -2{u - Ud)Vud - pVf, (3.26) 

51 = —Vp (g) MVu — Vu (g) M^^Vp + {MVu • Vp + up — fp + {u — Ud)^)!- (3.21) 


In the simple ease where M = I, assuming u,p ^ 67^(12), we know from the previous diseussion that (3.3) 
is satisfied. Noting that 


div(Vp (g) Vu) = AuVp + {D‘^p)'^Vu, 
div(Vu (g) Vp) = ApVu + {D‘^u)'^Vp, 
V {Vu ■ Vp) = D^uVp + D^pVu, 


the relation 


— div(Si) + So = 0 in 12 


(3.28) 


is equivalent to 

(—Am + u — f)Vp + (—Ap + p + 2 (m — Ud))Vu = 0 in 12. 

Therefore, we observe that the fundamental relation (3.3) between the tensors Si and So eorresponds to the 
strong solvability of the state and adjoint state equation. 


4 Descent directions 

In this paper we are interested in numerieal methods for shape optimization problems of the type 

minJ(12), (4.1) 

oeT 


where ip C V{D) is the admissible set. Assume J : ip —R is shape differentiable at 12 C 12 C R'^. 
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Definition 4.1 (descent direction). The vector field 9 G Cc'^{D^ R'^) is called a descent direction for J at 17 
if there exists an e such that 

J($f(17)) < J(17) for all t G (0,e). 

If the Eulerian semiderivative of J at 17 in direction 6 exists and if it is a descent direction then by definition 

dj{n){e) < 0. (4.2) 

Descent directions are used in iterative methods for finding approximafe (possibly local) minimizers of 
J(17). Typically, af a given sfarfing poinf 17, one defermines a descenf direction 9 and proceeds along fhis 
direcfion as long as fhe cosf funcfional J reduces sufficienfly using a sfep size sfrafegy. In fhis secfion we give 
a general setting for compufing descenf direcfions in fhe framework of gradienf mefhods using fhe domain and 
boundary represenfafions of fhe shape derivafive according fo Theorem 1.4. We show how a descenf direcfion 
9 wifh any regularify s > 1 can be obfained by solving an appropriafe partial differenfial equafion. We 
also show how fo deal wifh bound consfrainfs on 0. In order fo develop a setting allowing fo define general 
descenf direcfions, we recall sufficienl condifions for fhe solvabilify of fhe following operator equafion 

A9 = f, 

where A : E —)> E' is an operafor befween a Banach space E and ifs dual E'. Sufficienl condifions for fhe 
bijecfivily of A are given by fhe Iheorem of Minly-Browder [39, p.364. Theorem 10.49]. 

Theorem 4.2 (Minly-Browder). Lei (E; || • ||e) be a reflexive separable Banach space and A ; E —)• E' a 
bounded, hemi-conlinuous, monofone and coercive operafor. Then A is surjective, i.e. for each / G E' Ihere 
exisfs 0 G E such fhal A9 = f. Moreover if A is sfriclly monofone Ihen if is bijecfive. 

Lei A : E —)> E' be an operator on a reflexive, separable Banach space E safisfying fhe assumptions of 
Theorem 4.2 wifh A(O)0 > 0 for n G E. Assume dJ(Q) can be extended to E^ if necessary; for simplicily we 
keep fhe same nolalion for fhe extension. Infroduce fhe bilinear form 

^:ExE^R, B{9,C) ■■={A9, Oe',e- (4-3) 

Consider fhe variafional problem: 

(VP) Lind 01 G E such fhal ^(01, C) = -dJ{n){C) for all C G E, (4.4) 

Then fhe solution 0i of (VP) is a descenf direcfion since dJ{Q,){9i) = —B{9i, 0i) < 0. 

In cerlain siluafions if is desirable fo have bound consfrainfs on fhe shape perfurbalion. This may be 
handled by considering fhe more general case of a variafional inequalily. Given a subsel iT C E wifh 0 G iT, 
consider fhe variafional inequalily: 

(VI) Lind 02 G iT such fhal 73(02 ,02 - C) < dJ{n){C - 02) for all ( G K. 

The solulion 02 of (VI) yields a descenf direction for J at Q since faking ^ = 0 G iT we gel 


dJ(f7)(02) < -73(02,02) <0. 

In view of Theorem 1.4, we choose E C H^{D) where s is such fhal dJ{Q) : —)■ R'^ is 

conlinuous. When E is a Hilbert space, one may idenfify E' wifh E. Therefore if B is bilinear, coercive, and 
conlinuous, Ihen Lax Milgram’s lemma ensures fhal (VP) has a unique solulion. Lor all ofher cases we may 
have fo use Theorem 4.2 or similar resulls. 
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Remark 4.3. (a) Let E := Hq{D, R'^), B{9, C) := D9 : D( dx and (s D. Then (4.4) reads: find 

9 G Hq{D, R'^), such that B{9, () = —dJ(fl~^)(C) for all ( G Hq{D, R'^). Under the assumption that 
OU'*' G C^, 9\q+ G and £ H^{D \ U+), Proposition 3.3 yields 

[ g C ■ nds = dJ{il.'^){() for all C G Rq(L>, R'^), where g =—{D9^n — D9~n) ■ n. 

J dU+ 

This shows that the restriction to 90+ of the obtained descent direction 9 is more regular than the 
function g. 

(b) Let n be an extension of the unit normal to 0+ in D. If 9 defined on is a descent direction then 
{9 ■ n)n|gQ+ is also a descent direction, for the tangential part of 9 does not influence the derivative. 
Indeed define 9r ■= 9 — {9 ■ n)n, then by Nagumo’s theorem (iJ(O+)(0T-) = 0 and thus dJ(O+)(0) = 
dJ(O+)((0 • n)n). However, 9 and {9 ■ n)n lead to different transformations of the domains in general, 
indeed the tangential term actually has an influence for large deformations, which means <I)f(U+) / 
^{e influence appears for instance in the shape Hessian. 

5 Electrical impedance tomography 

We consider an application of the results above to a typical and important interface problem: the inverse 
problem of electrical impedance tomography (LIT) also known as the inverse conductivity or Calderon’s 
problem [6] in the mathematical literature. It is an active field of research with an extensive literature; for 
further details we point the reader toward the survey papers [5, 10] as well as [31] and the references therein. 
We consider the particular case where the objective is to reconstruct a piecewise constant conductivity a 
which amounts to determine an interface r+ between some inclusions and the background. We refer the 
reader to [2, 7,8,11,22,25,29] for more details on this approach. 

The main interest of studying LIT is to apply the approach developed in this paper to a problem which 
epitomizes general interface problems and simultaneously covers the entire spectrum of difficulties encoun¬ 
tered with severely ill-posed inverse problem. 


5.1 Problem statement 


Let D C R'^ be a Lipschitz domain, and U (Z D open sets such that D = U+ U U U r+, where 
r+ = 9U+ = u+ n fF and r = 9D = dn- \ r+; see Figure 1. In this section n denotes either the 
outward unit normal vector to D or the outward unit normal vector to U+. Decompose T as T = F u r„. 
Let a = cr+xn+ + where are scalars and / = /+Xn+ + f~Xn- where G H^{D). Consider 

the following problems: find Un G H\{D) such that 


/ aVun ■ Vz = fz + j gz for all z G H\{D) 

Jd Jd Jr„ 


and find Ud G such that 


[ aVud ■Vz= [ fz for all z G H^{D) 

Jd Jd 


(5.1) 


(5.2) 


H^D) 

hL{d) 

H^iD) 


{u G H^{D) I u = 0 on T^}, 

{v G H^{D) I u = 0 on T^, v = h on T^}, 
{v G H^{D) I u = OonT} 


where 
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and g G H represents the input, in this case the electric current applied on the boundary and h G 

is the measurement of the potential on r„, or the other way around, i.e. h can be the input and g 
the measurement. Define also the space 

PH^{D) := {u = u+xn+ + u~Xn -1 e n" G 

Consider the following assumption which will be used only to derive the boundary expression of the shape 
derivative but is not required for the domain expression: 

Assumption 5.1. The domains are of class C^, f G g G ^ (D) and 

h G for A: > 2. 

Applying Green’s formula under Assumption 5.1, equations (5.1) and (5.2) are equivalent to the following 


transmission problems where Un = u^Xn+ + Xn- 

and Ud = xq+ + Xn- ■ 


-cr+An+ = / in D+, 

-a~Au~ = fin n~, 

(5.3) 

u~ = OonTrf, 


(5.4) 

a~dnu~ = gon r„, 


(5.5) 

—(T^Au^ = / in D’*’, 

—a~Au'2 = / in 

(5.6) 

= OonTrf, 


(5.7) 

= hon Tn, 


(5.8) 


with the transmission conditions 

CJ+On'Wn = Cr~dnU~, = (J~dnU^ On T^, 

-I- — -\- — T^~l” \ * / 

Un=Un, OnT^. 

On Trf we impose homogeneous Dirichlet conditions, meaning that the voltage is fixed and no measurement 
is performed. One may take = 0 , in which case (5.1) becomes a pure Neumann problem and additional 
care must be taken for the uniqueness and existence of a solution. The situation T^ / 0 corresponds to partial 
measurements. Alternatively, it is also possible to consider a slightly different problem where each function 
Un and Ud has both the boundary conditions (5.5) and (5.8) on different parts of the boundary. 

Several measurements can be made by choosing sets of functions {gi}l^i and Writing Un,i and 

Ud,i for the corresponding states, the problem of electrical impedance tomography is: 

(EIT): Given {gi}l^i and {hi}l^i, find a such that Un,i = Ud,i in D ^or i = 1,.., I. 


(5.10) 
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Note that Un,i = Un,i{^~^) and Ud,i = actually depend on fl"'' through a = cr(n''“), however 

we often write Un,i and Ud^i for simplicity. In this section, we assume that the conductivities a~) are 
known, therefore the EIT problem (5.10) reduces to the following shape optimization problem where fl'*’ is 
the unknown. 

Given {fifJLi, and (o-+,cj-) with a = (J+xn+ + cr~Xn-^ 

find ri"'' such that Un,i = Ud^i in 17 for z = 1,.., I. 


Actually, the result for several measurements can be straightforwardly deduced from the case of one mea¬ 
surement by summing the cost functionals corresponding to each measurement, therefore in this section we 
stick to the case I = 1 of one measurement g for simplicity of presentation. In section 7 we consider several 
measurements for the numerics. 

The notion of well-posedness due to Hadamard requires the existence and uniqueness of a solution and the 
continuity of the inverse mapping. The severe ill-posedness of EIT is well-known: uniqueness and continuity 
of the inverse mapping depend on the regularity of a, the latter being responsible for the instability of the 
reconstruction process. Additionally, partial measurements often encountered in practice render the inverse 
problem even more ill-posed. We refer to the reviews [5, 10,31] and the references therein for more details. 
A standard cure against the ill-posedness is to regularize the inverse mapping. In this paper the regularization 
is achieved by considering smooth perturbations of the domains 17+. 

To solve problem (5.11), we use an optimization approach by minimizing the shape functionals 

= 11- Un{n+))^, (5.12) 

J2(f7+) = ^ - hf. (5.13) 

Since Ud,Un G H^{D) and h G 77^/^ (r„), Ji and J 2 are well-defined. Nofe fhaf Ji and J 2 are redundanf 
for fhe purpose of fhe reconstrucfion buf our aim is fo provide an efficienf way of computing fhe shape 
derivative of fwo functions which are offen encounfered in fhe liferafure. To compufe fhese derivatives we use 
fhe approach described in Secfion 2. Eirsf of all infroduce 

Fl{^d,^n) ■= \ I (5.14) 

^ JD 

F2{<^n) ■=\ [ - hf. (5.15) 

J 

Note fhaf Ji(17+) = rt„(n+)) and J 2 (n+) = F 2 {un{^'^))- Nexf consider ^ a subsef of V(D) 

and fhe Eagrangian C : x Hj{D) x H}{D) x H^{D) x Hj{D) R: 


£(Q+,¥?,'0) := alFl{^pd, ifn) + a2F2{^Pn) 


+ 

+ 


/ aV(fd-^'tpd- f'ipd+ / -O' dn't/Jdi^d - h) 

Jd JFn 


Id 


aVifn ■ VV'n - fi’n “ / 

7r„ 


(5.16) 


where ip := {ipd,ipn) and ij: := {^di'i’n)- The term —a dnipd in the second integral of (5.16) is used fo 
enforce fhe boundary condifion (5.8). Infroduce fhe objective funcfional 


J(17+) := ai Ji(f7+) -f a2J2(^7+). 


To compufe fhe derivafive of fhe Eagrangian depending on (5.3)-(5.8) we apply fhe averaged adjoinf 
mefhod from Secfion 2. 
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5.2 State and adjoint equations 

The state u := {ud,Un) and adjoint state p := {pd,Pn) are solutions of the equations: 

u, p)(V') = 0 for all ■0 G Hq{D) x H\{D), 
u, p)((^) = 0 for all p G H\{D) x H\{D). 


(5.17) 

(5.18) 


Writing (5.17) explicitely, one ean obtain easily the state equations (5.1) and (5.2). Then (5.18) yields the 
equation for the adjoint pd'- 


5^^£(n+, u, p){pd) = 0, for all 0d G 


whieh leads to 


Id 


aVpd-Vpddx = -ai / {ud - Un)pddx - / -a dnPdPdds for all G 77](77) (5.19) 

Jd 7r„ 


whieh is the variational formulation for the adjoint state pd- This yields the following variational formulation 
when test funetions are restrieted to 77] (77): 


/ aVpd ■ Vip dx =—ai / {ud — Un)(pdx for all G 77] (77). 
Jd Jd 


(5.20) 


Ifweuse Assumption 5.1, we getprf G PH^{D) and using Green’s formula in 17+ and 17 with(/9 G C'“(17+) 
and p G (17“), we obtain the strong form 


— dw{aVpd) = —ai{ud — Un) in 17+ and 17 . 


(5.21) 


Henee using now Green’s formula in (5.19) and using (5.21) gives 

/ [(ydnPd]v+‘^d ds+ {adnPd - (y~dnPd)Pd ds = d for all (pd G 77] (77), 

7r+ 7r„ 

where [(TdnPd\T+ = o'~^dnPd — o'~dnPd is the jump of adnPd aeross r+. Sinee the integral on T^ above 
vanishes and pd G Hq{D), we obtain 

Pd = 0on r, (5.22) 

<y"^dnPd = cr~dnPj on r+. (5.23) 

Finally solving 

9^„£(17+,u,p)((^„) = 0, for all G H^{D), 
leads to the variational formulation 



Un)Pn + Cr'^Pn ' V(^n + 



h)(pn = 0 


for all ifn G 77] (77). 

Similarly as for we get, under Assumption 5.1, p„ G PH^{D) and the strong form 

— div(crVpn) = C(i{ud — Un) in 17+ and 17“, 
crdnPn = -a2{un - /i) on T^, 

Pn = Oon r^i, 

<y^dnpt = o-~dnPn on 7+, p+ = p~ on r+. 


(5.24) 


(5.25) 

(5.26) 

(5.27) 

(5.28) 
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5.3 Shape derivatives 

Let us consider a transformation <l>f defined by (1.3) with 9 G Cl{D,Y{f‘). Note that ^t{D) = D but in 
general <l>f (fl"'') / We use the notation 11+ {t) := (11+). Our aim is to show the shape differentiability 

of J(ll+) with the help of Theorem 2.1. Following the methodology described in Section 2 we introduce 

:=£(0+(t),(^oT>-\V’o^r^)- (5-29) 

We proceed to the change of variables = y in (5.29) to get the canonical form (2.1). First of all let 

us denote /f^+(t) = /+Xn+(t) + f~Xn-(t) and cjf^+p) = a+xn+{t) + cr~Xn-{ty^ recall that f 7 + are scalars 
but /+ are functions. Then note that the change of variables <l>t(x) = y leads to considering the following 
functions inside the integrals: 

o-Q+p) o<^t = o-^Xn+it) o^t + (y~Xn-(t) °^t = cr^Xn+ + cr~Xn- = 

fn+{t) o^t = Xn+{t) o^t + r o ^tXn-{t) o o ^tXn+ + f~ ° ^tXn-- 

Thus we introduce the function ft := /+ o ^t Xn+ + f~ ° Xn- ■ Now we obtain the canonical form (2.1) 
for the Lagrangian: 


= a{t,if,'il;) + h{t,^p), (5.30) 

with 

a{t, ip, V’) := / aA{t)Vpd ■ - Mdf.it) - / a~^dni’d{Vd - h) 

JD JVn 

+ / aA{t)Vpn ■ VV’n - fti’nf.it) “ / gf’n-, 

JD JFn 

b{i,V>)- = Y [ iv>d- PnfCit) + Y [ iv>n-hf 

J D J ji 

where the Jacobian f^{t) and A{f) are defined as f,{t) := det(i9<l>t) and A{t) := . In fhe 

previous expression (5.30), one should note fhaf fhe infegrals on subsefs of F are unchanged since = I 
on r. Thus we have <l>f {D) = D, however fhe terms inside fhe infegrals on D are modified by fhe change of 
variable since / I inside D. Note fhaf 

J(52+(t)) = G{t, u*, V^), for all V’ G hI{D) x hI{D), 

where u* = (rr^, n^) ;= {un,t ° Ud,t ° and Un,t, Ud,t solve (5.1),(5.2), respectively, wifh fhe domains 
12+ and 12“ replaced by 12+(f) and 12“ (f). As one can verify by applying a change of variables fo (5.1) and 
(5.2) on fhe domain 12+(t) fhe functions satisfy 


[ aA{t)Vui-Vi>n= [ Mn+ [ gi’n for all iJn G HdiD), 
JD JD Jr„ 

[ aA{t)Vu\-V^d= [ Md for allied eH^{D). 


(5.31) 

(5.32) 


JD JD 

Applying sfandards esfimafes for elliptic parfial differential equafions and fhe facl fhaf A{t) is uniformly 
bounded from below and above for t small enough, we infer from equations (5.31),(5.32) fhe existence of 
consfanfs Ci, 6*2 > 0 independenf of t and r > 0 such fhaf for all t G [ 0 ,r]: 

\\'>4i\\h^{d) < Gi, and \\un\\m(D) < C' 2 - 


(5.33) 
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From these estimates, we get ^ Wd and u\^ Wn in H^{D) as t —)■ 0. Passing to the limit in (5.31) and 
(5.32) yields Wd = Ud and Wn = Un by uniqueness. 

Let us now eheek Assumption (HO) and the eonditions of Theorem 2.1 for the funetion G given by (5.30) 
and the Banaeh spaees E = H\{D) x H^{D) and F = Hq{D) x H^{D). First of all equation (2.2) admits 
a unique solution u* ;= (n^, «^) for eaeh t G [0, r]. The eonditions of Assumption (HO) are readily satisfied 
and also the funetion G is affine with respeet to V’ = (V'd, ipn)- 

Regarding the eonditions of Theorem 2.1, first note that applying Lax-Milgram’s lemma, we eheek that 
both equations (5.34) and (5.35) have indeed a unique solution in F = H^{D) x HjiD): 

[ aA{t)Vpd-'^‘Pd + ^ [ ^{t){Ud + Ud- {ui + Un))pd- [ (T~^dnPd^d = 0, (5.34) 

Jd ^ Jd Jr„ 

[ aA{t)Vpi-V(pn-^ [ ^{t)iUd + Ud- {ui + Un))(pn + ^ [ + Un - 2h)(pn = 0, (5.35) 

Jd ^ Jd ^ Jv„ 

for all (pd, ‘Pn in H^{D). Therefore there exists a unique solution p* = {ph,Pd) averaged adjoint 

equation (2.3). 

Now we eheek Assumption (HI). Testing (5.34) with pd = Pd nnd (5.35) with pn = Pn, we eonelude by 
an applieation of Holder’s inequality together with (5.33) the existenee of eonstants Gi, G 2 and r > 0 sueh 
that for all t G [0, r] 

\\Pd\\H^{D) < C*!) and \\Pn\\m{D) < C 2 . 

We get that for eaeh sequenee eonverging to zero, there exists a subsequenee also denoted sueh that 
Pd" Id and pI^ qn for two elements Qd, qn £ Passing to the limit in (5.34) and (5.35) yields 

Qd = Pd and qn = Pn by uniqueness, where pd and pn are solutions of the adjoint equations. Sinee the limit 
is unique, we have in faet ^ pd and p\^ ^ pn^^t ^ 0. Finally, differentiating G with respeet to t yields 

dtG{t,p,'ip) = Y j^i‘Pd-‘fnfC{t)tT{D9tD^^^) 

+ [ aA'{t)Vpd • VV’d - Mdm - i^dVft • etm 

Jd 

+ [ aA'{t)Vpn • VlPn - Mnm tr(F0tF4>,-l) - V'nV/i • Otm- 

Jd 


where 

^ft ■■= V/+ o ^tXn+ + V/“ o ^tXn-, 9t = 0 o 

Aft) = tT{D9^D<^^fA{t) - D<^^^D9tA{t) - D9tA{t))^ 

and D9t is the Jaeobian matrix of 9t. In view of 0 G (7^(1), R'^), the funetions t i-G D9t and t 1 —)■ 
tr(F)0t4>jr^) = div(0) o are eontinuous on [0, T]. Moreover pd, Pn, fdi fn are in H^{D), f G PH^{D) 
so that dtG{t,p,'ijj) is well-defined for all t G [0,r]. Using the weak eonvergenee of p* and the strong 
differentiability oft ^ A{t) and 1 1 —)■ ^{t) it follows 


^.^^ G(f,uO,p*)-G(0,uO,p*) 
t\o t 


5tG(0,u°,p°). 


(5.36) 


Thus we have verified all assumptions from Theorem 2.1. This yields 

dJ(U+)(0) = ^ (G(f, u*, tf)) |t=o = dtGf), uO, pO) for all V’ G F = hI{D) x Fi(F), 


and therefore we have proved the following result. 
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Proposition 5.2 (distributed shape derivative). Let D C R'^ be a Lipschitz domain, 9 G Cl{D,R% 
f G PH\D), g G R-i/2(r„), h G 52+ C is an open set, then the shape derivative of 

J(52+) is given by 


dJ(52+)(6') = j - '“n)^ - f{Pn+Pd)j div^ 

+ [ -iPd+Pn)'^f ■9 + aA'{0){Vud-Vpd + Vun-Vpn), 


(5.37) 


ID 


where V/ := V/+ Xq+ + V/ Xn-> ^^(0) = (div6l)/ — 720^ — D9, Un, Ud are solutions of (5.1),(5.2) and 
Pn,Pd of (5.24), (5.19). 

The shape derivative (5.37) also has the tensor representation corresponding to (3.1): 


where 


(iJ(52+)(0) = [ Si-09 +So-9, 
Jd 


(5.38) 


Si = -a{Vud <8) Vpd + Vpd ® Vud + Vu^ (g) Vpn + Vpn ® Vtt„) 

+ a{VUd ■ S/pd + S/Un ■ S/pn)I + “ /(Pn +Pd)) 

So = -(Pd +Pn)v/. 

Note that the volume expression of the shape gradient in Proposition 5.2 has been obtained without any 
regularity assumption on 52+. In order to obtain a boundary expression on the interface r+ we need more 
regularity of 57+ provided by Assumption 5.1. If it is satisfied, we can apply Corollary 3.5 to obtain directly 
the boundary expression of the shape derivative, using mainly the standard tensor relation (Vud <S> Vpd)n = 
(Vpd ■ n)'Vud, which yields Proposition 5.3. 

Proposition 5.3 (boundary expression). Under Assumption 5.1 and 9 G (72, R*^) the shape derivative of 
J(52+) is given by 

d<7(57 )(9) — / [it( djiUddfiPd dnUn9nPn)]Y'+ 9 • n 
Jt+ 

+ / ([o']r+(Vr+tid • Vr+Pd + Vr+Un • Vr+Pn) - [/]r+(Pn+Pd))6'• n. 

7r+ 

Note that our results cover and generalize several results that can be found in the literature of shape 
optimization approaches for EIT, including [2,22]. For instance taking a 2 = = 0 in Proposition 5.3 we 

get pd = 0 which yields the same formula as the one obtained in [2, pp. 533]. 

Note also that from a numerical point of view, the boundary expression in Proposition 5.3 is delicate 
to compute compared to the domain expression in Proposition 5.2 for which the gradients of the state and 
adjoint states can be straightforwardly computed at grid points when using the finite element method for 
instance. The boundary expression, on the other hand, needs here the computation of the normal vector and 
the interpolation of the gradients on the interface r+ which requires a precise description of the boundary 
and introduces an additional error. 


6 Level set method 

The level set method, originally introduced in [35], gives a general framework for the computation of evolving 
interfaces using an implicit representation of these interfaces. The core idea of this method is to represent the 
boundary of the moving domain 52+(f) C 72 G R^ as the level set of a continuous function i;()(-, f) : 72 —)> R. 


22 


Let us consider the family of domains C D as defined in (1.4). Each domain can be defined 

as 

:= {x G 22, (j){x,t) < 0} (6.1) 

where cj) : D x R+ ^ R is continuous and called level set function. Indeed, if we assume | V(/)(-, f)| / 0 on 
the set {x G D, (j){x, t) = 0} then we have 

= {x £ D, 4>{x,t) = 0}, (6.2) 

i.e. the boundary f)n'''(f) is the zero level set of </>(-, f). 

Let x{t) be the position of a particle on the boundary f)n+(f) moving with velocity x{t) = 6{x{t)) 
according to (1.3). Differentiating the relation (/)(x(f),f) = 0 with respect to t yields the Hamilton-Jacobi 
equation: 

dt4>{x{t),t) + 6{x{t)) ■ \7(l){x{t),t) = 0 in x R'*", 

which is then extended to all of D via the equation 

dtf{x,f) + 9{x) ■V(l){x,f) = ini2xR^, (6.3) 

or alternatively to ?7 x R"'' where ?7 is a neighbourhood of 

Traditionally, the level set method has been designed to track smooth interfaces moving along the normal 
direction to the boundary. Theoretically, this is supported by Theorem 1.4, i.e. if the domain and the 

shape gradient are smooth enough then the shape derivative only depends on 0 • n on (9D+(f). In this case, 
we may choose for the optimization a vector field 9 = iDniT- on 012+(f). Then, noting that an extension to D 
of the unit outward normal vector n to 12+(t) is given byn = V(/)/|V(/>|, and extending to all of D, one 
obtains from (6.3) the level set equation 

0t(/) + t2n|V(/>| = 0 inZ2xR+. (6.4) 


The initial data (/)(x, 0) = fo{x) accompanying the Hamilton-Jacobi equation (6.3) or (6.4) is chosen as the 
signed distance function to the initial boundary 012+(0) in order to satisfy the condition |Vri| / 0 on 012+, 
i.e. 


f d(x, 012+(0)), if X G (!2+(0))'^, 
\ —d(x, 012+(0)), ifxGl2+(0). 


(6.5) 


6.1 Level set method and domain expression 

In the case of the distributed shape derivative, for instance (5.37) or (5.38), f is not governed by (6.4) but 
rather by the Hamilton-Jacobi equation (6.3). Indeed we obtain a descent direction 9 defined in D by solving 
(4.4), where (2J(12+) is given by Proposition 5.2 which can subsequently be used in (6.3) to compute the 
evolution of f. On the other hand, in the usual level set method, one solves a PDE on the boundary 912+ 
in an analogous way as for (4.4) (for instance using a Laplace-Beltrami operator), and uses the boundary 
expression from Proposition 5.3 to obtain dn = 9 • n on 912+. 

Numerically it is actually more straightforward in many cases to use (6.3) instead of (6.4). Indeed, when 
using (6.4), dn is initially only given on 912+(2) and must be extended to the entire domain D or at least to a 
narrow band around 912+(2). Therefore it is convenient to use (6.3) with 9 already defined in D as is fhe case 
of the distributed shape derivative, which provides an extension to Z2 or to a narrow band around 912+(2). 

In shape optimization, usually depends on the solution of one or several PDEs and their gradient. 
Since the boundary 912+(2) in general does not match the grid nodes where f and the solutions of the partial 
differential equations are defined in the numerical application, the computation of requires the interpola¬ 
tion on 912+(2) of functions defined at the grid points only, complicating the numerical implementation and 
introducing an additional interpolation error. This is an issue in particular for interface problems where is 
the jump of a function across the interface, as in Proposition 5.3, which requires multiple interpolations and 
is error-prone. In the distributed shape derivative framework 9 only needs to be defined at grid nodes. 



23 


6.2 Discretization of the level set equation 

Let D be the unit square D = (0,1) x (0,1) to fix ideas. For the discretization of the Hamilton-Jacobi equation 
(6.3), we first define the mesh grid corresponding to D. We introduce the nodes Pij whose coordinates are 
given by {iAx,jAy), I < i,j < N where Ax and Ay are the steps discretization in the x and y directions 
respectively. Let us also write = kAt the discrete time for k £ N, where At is the time step. We are 
seeking for an approximation ~ 

In the usual level set method, the level set equation (6.4) is discretized using an explicit upwind scheme 
proposed by Osher and Sethian [34,35,40]. This scheme applies to the specific form (6.4) but is not suited to 
discretize (6.3) required for our application. Equation (6.3) is of the form 

dt(j) + = 0 in D X R+. (6.6) 


where H{V(j)) := 9 ■ V</> is the so-called Hamiltonian. We use the Local Lax-Friedrichs flux originally 
conceived in [36] and which reduces in our case to: 


H 


LLF/- 


{p ,P ,9 ,Q )= H 


p +p^ q +q^\ 1 . + ^ 1 . + y 

-^{P -P )« -i^KP -P 


where a* = \6x\, = \6y\, 6 = Oy) and 

4^ij 


P — Px 4’ij 
q~ = Dy4>ij 


Ax 

Ay 


P'^ = = 

q+ = D+ct)ij = 


4^' 




^ Ax 

4^1, j+i 4^ij 

Ay 


are the backward and forward approximations of the x-derivative and y-derivative of 4> at Pij, respectively. 
Using a forward Euler time discretization, the numerical scheme corresponding to (6.3) is 


4 ^' = 4 - - H^^^{P~,P'^, q~,q^) 




(6.7) 


For numerical accuracy, the solution of the level set equation (6.3) should not be too flat or too steep. 
This is fulfilled for instance if 4> is the distance function i.e. |V(^| = 1. Even if one initializes 4> using a 
signed distance function, the solution 4> of the level set equation (6.3) does not generally remain close to a 
distance function. We may occasionally perform a reinitialization of 4> by solving a parabolic equation up to 
the stationary state; see [16,17,38]. Although in the level set method this reinitialization is standard, in the 
case of the distributed shape gradient, we observe experimentally that the level set function 4> stays close to 
a distance function during the iterations and we do not need to reinitialize. The regularization of the shape 
gradient could explain this observed stability of the level set function. 

The computational efficiency of the level set method can be improved by using the so-called “narrow 
band” approach introduced in [1], which consists in computing and updating the level set function only on a 
thin region around the interface. This allows to reduce the complexity of the problem to N log{N) instead of 
N'^ in two dimensions. In this paper we do not implement this approach but we mention that it could also be 
applied to the distributed shape derivative approach and equation (6.3) by taking 6 with a support in a narrow 
band around the moving interface, which can be achieved by choosing the appropriate space E in (4.3). 


7 Application and numerical results 

7.1 Electrical impedance tomography 

In this section we give numerical results for the problem of electrical impedance tomography presented in 
Section 5.1, precisely we look for an approximate solution of the shape optimization problem (5.11). Using 
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the notations of Section 5.1 we take D = (0,1) x (0,1) and = 0, i.e. we have measurements on the entire 
boundary F. For easiness of implementation, we consider a slightly different problem than the one in Section 
5.1. Denote Ft, F;,, F/ and F^ the four sides of the square, where the indices t, b, I, r stands for top, bottom, 
left and right, respectively. We consider the following problems: find G 


/ aVun-V(p= fif+ gif for all (/? G Ffotb(^) 
Id Jd JriuTr 


and find Ud G such that 


'D 


where 


aVud • 

S/ip = 

■f 

fV’+ [ 

9^ 

for all 

f^K 



Jd 


^rtUF;, 



K 

(D) 

= {r; 

G 

H\D)\ 

V = 

on Ft 

UFb}, 

Hi 

(D) 

= {r; 

G 

H\D)\ 

V = 

/i on F; 

UFJ, 

Hltb 

(D) 

= {n 

G 

H\D)\ 

V = 

0 on Ft 

UFfe}, 

Kir 

(D) 

= {n 

G 

H\D)\ 

V = 

0 on F; 

UFJ. 


(7.1) 


(7.2) 


In our experiments we choose / = 0. The results of Section 5.1 can be straightforwardly adapted to equations 
(7.1), (7.2). 

We use the software package FEniCS for the implementation; see [30]. The domain D is meshed using 
a regular grid of 128 x 128 elements and we describe the evolution of the interface F+ using the level set 
method from Section 6. The conductivity values are set to uo = 1 and ai = 10. 

We obtain measurements hi corresponding to fluxes gi,i = 1,by taking the trace on F of the solution 
of a Neumann problem where the fluxes are equal to pj. To simulate real noisy EIT data, the measurements 
hi are corrupted by adding a normal Gaussian noise with mean zero and standard deviation J * ||/ii||oo 5 where 
(5 is a parameter. The noise level is computed as 


noise = 


J2i=i \\hi - hi\\L2(^Y') 


(7.3) 


where hi is the noisy measurement and hi the synthetic measurement without noise on F. 
We use a variation of the functional (5.12), i.e. in our context: 


j(D+) = V/Ti [ l{ud,iin+)-un,i{n+))^, 
i=i ^ 


(7.4) 


where Ud,i and Un,i correspond to the different fluxes gi. Here the coefficients fi are weights associated to 
the fluxes gi. In our experiments we choose the weights m such that each term of the sums in (7.4) are 
equal to 1 on initialization in order to have a well-distributed influence of each term. Practically, the /r* are 
thus calculated during the first iteration. We use the distributed shape derivative dJ(D+) from Proposition 
(5.2). We obtain a descent direction by solving (4.4) with E a finite dimensional subspace of Hq{D) and 
B{v, w) = fjy Dv ■ Dw. We choose E to be the space of linear Eagrange elements. 

Since we use a gradient-based method we implement an Armijo line search to adjust the time-stepping. 
The algorithm is stopped when the decrease of the functional becomes insignificant, practically when the 
following stopping criterion is repeatedly satisfied: 

J(D+) - J(F!++,) < 7(J(F!+) - J(F!+)) 
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Figure 2: Reconstruction (continuous contours) of two ellipses (dashed contours) with different noise levels 
and using three measurements. From left to right and top to bottom: initialization (continuous contours - top 
left), 0% noise (367 iterations), 0.43% noise (338 iterations), 1.44% noise (334 iterations), 2.83% noise (310 
iterations), 7% noise (356 iterations). 


where denotes the /c-th iterate of We take 7 = 5.10“® in our tests. 

In Figure 2 we compare the reconstruction for different noise levels computed using (7.3). We take in this 
example 1 = 3, i.e. we use three fluxes gi, i = 1,2, 3, defined as follows: 

Pi = 1 on F/ U Fr and pi = —1 on F^ U F;,, 

P 2 = 1 on F; U Fi and p 2 = —1 on F^ U F;,, 

P 3 = 1 on Fz U Ffo and gs = -1 on F,, U F*. 

Without noise, the reconstruction is very close to the true object and degrades as the measurements become 
increasingly noisy, as is usually the case in FIT. However, the reconstruction is quite robust with respect 
to noise considering that the problem is severely ill-posed. We reconstruct two ellipses and initialize with 
two balls placed at the wrong location. The average number of iterations until convergence is around 340 
iterations. 

In Figure 3 we reconstruct three inclusions this time using 1 = 7 different measurements, with 1.55% 
noise. The reconstruction is close to the true inclusion and is a bit degraded due to the noise. Figure 4 shows 
the convergence history of the cost functional in log scale for this example. 

Our algorithm gives good results in comparison to existing results in the literature using level set methods 
to solve the BIT problem. In [II] the BIT problem has been treated numerically using a level set method, 
which is not based on the use of shape derivatives but on the differentiation of a smooth approximation of the 
Heaviside function to represent domains. In [22] the level set method using the boundary expression of the 










26 


shape derivative is used based on equation (6.4). 

Our algorithm converges fast in comparison to [11,22]: convergence occurs after around 300 iterations. 
In [11] convergence occurs between 200 iterations for one inclusion and up to 50000 iterations for two 
inclusions. In [22] convergence occurs after 2000 or 10000 iterations on two examples with three inclusions. 
Concerning measurements we obtain good reconstruction of two inclusions with 1 = 3 and three inclusions 
with 1 = 7 , while in [11] sets of 4,12, 28 and 60 measurements are used but usually 60 measurements are 
required for complicated shapes such as two inclusions. In [22], 60 measurements are used. Nevertheless, 
our results are not directly comparable since the conductivities are unknown in [11,22], which makes the 
inverse problem harder and might explain the slower convergence. Also, the reconstructed shapes are not the 
same although the complexity of the unknown shapes is comparable since we also consider two and three 
inclusions as in [11,22]. Only an exact comparison using the same problem, test case, initialization, noise 
level, number and type of measurements could allow to conclude. 



Figure 3: Initialization (continuous contours - left) and reconstruction (continuous contours - right) of two 
ellipses and a ball (dashed contours) with 1.55% noise (371 iterations) and using seven measurements. 



Figure 4: History of cost functional corresponding to Figure 3 in logarithmic scale. 
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